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Molecular hydrogen emission from protoplanetary disks 
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Abstract. We have modeled self-consistently the density and temperature profiles of gas and dust in protoplanetary 
disks, taking into account irradiation from a central star. Making use of this physical structure, we have calculated 
the level populations of molecular hydrogen and the line emission from the disks. As a result, we can reproduce 
the observed strong line spectra of molecular hydrogen from protoplanetary disks, both in the ultraviolet (UV) 
and the near-infrared, but only if the central star has a strong UV excess radiation. 
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1. Introduction 

It has long been believed that infrared excesses over the 
stellar photospheric emissions, often observed in the spec- 
tral energy distributions of young stellar objects (YSOs), 
arise from dusty circumstellar disks (Cohen 1974; Kenyon 
& Hartmann 1987; Beckwith & Sargent 1993). More di- 
rect evidence for circumstellar disks has been found re- 
cently in the form of optically thick dust lanes against the 
scattered light of the surrounding optically thin nebulae 

. at optical and near- infrared wavelengths (e.g., Burrows et 
al. 1996; Stapelfeldt et al. 1998, 2003; Cotera et al. 2001; 
Jayawardhana et al. 2002) and as a near-infrared image of 

1 the scattered light of the disk itself (e.g., Fukagawa et al. 
2004). 

Furthermore, thanks to recent high spectral resolution 
and high sensitivity observations, it has become possible 
to detect various molecular line emission from protoplan- 
etary disks (e.g., Dutrey et al. 1997). In particular, obser- 
vation of molecular hydrogen line emission from the disks 
is important because it directly traces gaseous masses of 
the disks, which are connected with giant planet forma- 
tion, without assuming the dust-to-gas ratio or CO-to-H2 
ratio. Pure rotational molecular hydrogen line emission 
was detected towards YSOs and debris-disk objects with 
the Infrared Space Observatory (ISO) (Thi et al. 1999, 
2001a, 2001b), while no emission has yet been observed 
with ground-based telescopes (Richter et al. 2002; Sheret 
et al. 2003; Sako et al. 2005). Ro- vibrational molecular hy- 
drogen line emission with narrow line widths was detected 
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towards some T Tauri stars (Weintraub et al. 2000; Bary 
et al. 2002, 2003; Itoh et al. 2003). In addition, fluores- 
cent molecular hydrogen line emission in the ultraviolet 
(UV) wavelength band was observed towards some classi- 
cal T Tauri stars (Herczeg et al. 2002, 2004; Bergin et al. 
2004). However, at present there is no established model 
for this emission which takes into account the global phys- 
ical structure of protoplanetary disks. 

Historically, molecular hydrogen emission has been ob- 
served towards various kinds of astronomical objects, such 
as shock surfaces associated with star forming regions, 
reflection nebulae illuminated by nearby massive stars, 
planetary nebulae, supernova remnants, external galax- 
ies, etc. (e.g., Beckwith et al. 1978; Brown et al. 1983; 
Hasegawa et al. 1987; Burton et al. 1992), and studied the- 
oretically, for example, under the conditions of shock or 
photon-dominated regions (e.g., Black & Dalgalno 1976; 
Shull 1978; Hollenbach & McKee 1979; Draine et al. 1983; 
Pineau des Forets et al. 1986; Black & van Dishoeck 
1987; Wagenblast & Hartquist 1988; Sternberg 1988, 1989; 
Sternberg & Dalgalno 1989; Draine & Bertoldi 1996). 
These studies have suggested the importance of physical 
condition of the objects - density, temperature, and UV 
irradiation - for exciting molecular hydrogen. 

Now, the physical structure of protoplanetary disks 
is thought to be controlled by irradiation from the cen- 
tral star. Modelling the density and temperature pro- 
files of dust and gas in the disks has been developed in 
more and more realistic ways as more detailed observa- 
tional data become available (e.g., Kusaka et al. 1970; 
Kenyon & Hartmann 1987; Chiang & Goldreich 1997; 
D'Alessio et al. 1998; Nomura 2002; Dullemond et al. 2002; 
Dullemond & Dominik 2004; Kamp & van Zaldelhoff 2001; 
Gorti & Hollenbach 2004; Glassgold et al. 2004; Kamp & 
Dullemond 2004; Jonkheid et al. 2004). 
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In this paper, we have modeled the density and tem- 
perature profiles of protoplanetary disks self-consistently, 
taking into account the irradiation from the central T 
Tauri star. We have used these to investigate the abun- 
dance and excitation of molecular hydrogen in the disk, 
and the observational properties of the molecular hydro- 
gen emission from the disk. In the following section, we 
present the physical model of the disk:- the density and 
temperature profiles of the gas and dust in the disk on 
the assumptions of vertical hydrostatic equilibrium and 
local thermal and radiative equilibrium. In Sect. 3, we 
calculate the abundance and level populations of molecu- 
lar hydrogen, assuming statistical equilibrium among the 
levels. Making use of these physical and chemical profiles, 
we compute molecular hydrogen emission from the disk 
at infrared and ultraviolet wavelengths and compare with 
the observations in Sect. 4. Finally, the results are sum- 
marized in Sect. 5. 

2. Physical model 

We consider an axisymmetric disk surrounding a central 
star with the physical parameters of typical T Tauri stars; 
a mass of = 0.5M Q , a radius of i?* = 2R Q , and a 
temperature of T* = 4000K (e.g., Kenyon & Hartmann 
1995). 

2.1. Basic equations for the disk structure 

The gas temperature and density distributions of the disk 
are obtained self-consistently by iteratively solving the 
equations for hydrostatic equilibrium in the vertical di- 
rection and local thermal balance between heating and 
cooling of gas. The equation for vertical hydrostatic equi- 
librium is given in cylindrical coordinates (x, z) by 

dP 2 dp 

dz s dz 

where P and p represent the pressure and density, respec- 
tively. The sound speed c s is defined as c s 2 = dP/dp = 
kT/m^, where k and T represent Boltzmann's constant 
and the gas temperature, and the mean molecular mass 
m M is set to be m M = 2.3toh (toh is the hydrogen mass). 
The vertical gravitational force is set g z — GM*z/(x 2 + 
z 2 ) 3 / 2 , where G is the gravitational constant. The condi- 
tion, J z 2 °° p(x,z)dz — £(2), is set to be satisfied, where 
T,(x) is the surface density at a radius x defined below. 
We put p(x, Zoo) = 1.67 x 10~ 20 g cm~ 3 for the boundary 
condition. The equation for the detailed energy balance at 
each point in the disk is given by 



-P9z 



(1) 



r = A, 



(2) 



where T and A are the sum of the relevant gas heating and 
cooling rates, respectively (see Appendix A for details of 
the heating and cooling processes). 

The dust temperature profile plays an important role 
in determining the disk structure because the gas temper- 
ature is well coupled with the dust temperature near the 



midplanc of the disks where the density is high enough 
for efficient collisions between gas and dust particles (see 
Sect. 2.4). In this paper we obtain the dust temperature 
profile by assuming local radiative equilibrium between 
absorption and reemission of radiation by dust grains at 
each point in the disk in spherical coordinates (R, 0), 



4ir [ dv K „(R,<d)B u [T d (R,<3)} 
Jo 



-f 

Jo 



dvn l/ (R,Q) f dpd(j)I^(R,B; p, cf>), 



(3) 



where n u ,I v , B v (Td), and T d represent the monochromatic 
absorption coefficient, the specific intensity, the Planck 
function for blackbody radiation at a frequency v, and 
the dust temperature, respectively. Local thermodynamic 
equilibrium, t] v = n v B v (Td), is assumed, where r\ v is the 
monochromatic emissivity. The specific intensity is calcu- 
lated by solving the axisymmetric two-dimensional radia- 
tive transfer equation, 



T V (R, 9; /i, <j>) = [ k„(R', e')p(R', 9') 
Jo 



xB v [T d {R',Q')]e- T ^ R '> & 'Us', (4) 

where t v (R! , 9') is the specific optical depth from a point 
(R',0') to (R, 9), by means of the short characteristic 
method in spherical coordinates (Dullemond & Turolla 
2000). Here we neglect the effect of scattered light (cf. 
Dullemond & Natta 2003) . As heating sources we consider 
the radiative flux produced by the viscous dissipation (a- 
model) at the midplane of the disk, F v i s = (9/4)Y,ac s ^il, 
and the irradiation from the central star F^ar (Appendix 
C; see also Nomura 2002 for details of the model). 

The surface density distribution of the disks, S(x), is 
determined based on the standard accretion disk model 
(e.g., Lynden-Bell & Pringle 1974; Pringle 1981), by equat- 
ing the gravitational energy release of accreting mass to 
thermal heating via viscous dissipation at the disk mid- 
plane at each radius, x, 
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where c s0 and £1 = (GM r /x ) ' represent the sound 
speed at the midplane and the Keplarian frequency, re- 
spectively. Here, we adopt a = 0.01 for a viscous parame- 
ter, and assume that the disk has a constant mass accre- 
tion rate of M = 10~ 8 M Q yr -1 . The disk mass between 
the inner radius r; n = i?» = 2i? Q and the outer radius 
r out = 100AU is then 6 x 10" 3 M Q . 

2.2. Ultraviolet radiation fields 

Ultraviolet (UV) radiation fields affect the gas tempera- 
ture of the protoplanetary disks through the grain photo- 
electric heating (Appendix A.l) as well as the photodis- 
sociation and photoionization processes (Appendix B and 
Sect. 3). 
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Fig. 1. The gas temperature (a-1) (b-1) and density (a-2) (b-2) distributions in the x-z plane for the models (a) with 
and (b) without UV excess radiation of the central star. 



The main UV radiation sources are generally thought 
to be due to the central star and the interstellar radiation 
field. In this paper we simply adopt a 1+1 dimensional 
UV radiation field. In the radial direction, we calculate 
the radiation field which directly comes from the central 
star as 



F V . R {R, 0) = fF v , aiai exp(-T^ij), t v . r 



XvpdR, (6) 



R, 



where i^ jStar is the specific radiation field at the stellar 
surface, / = (7r/4)(i?„/i?) 2 the reduction of radiation by 
a geometrical effect. t v .r is the specific optical depth from 
the stellar surface (i?»,0) to a point (R, 0), and \v is 
the monochromatic extinction coefficient defined by the 
absorption and scattering (<t„) coefficients as \v = 
k v -\-<Jv I n the vertical direction, we consider the radiation 
source of the interstellar radiation field, i^isRF, plus the 
contribution of the scattering of radiation from the central 
star, 

F V , Z (X, z) = i^ISRF CXp[-T 1/jZ (z 00 )] 

+2ttuv / n v {x,z')p{x,z')F VtR {x,z')e~ T "^ z ' ) dz l 



f 

J z 



XvPdz", 



(7) 



where lo v is the monochromatic albedo defined by the ab- 
sorption and extinction coefficients as io v = ov/x„, and 
T v , z (z') is the specific optical depth from a point (x, z) to 
(x, z'). The approximate treatment of the scattered light 
in this model could overestimate the UV radiation fields 
in the disks. Fully two-dimensional radiative transfer cal- 
culation including the scattering process should be done 
in future (e.g., van Zadelhoff et al. 2003). 

Now, it is known observationally that many classical T 
Tauri stars have excess continuum radiation in the ultra- 
violet region, compared to the main-sequence stars of sim- 
ilar effective temperature (e.g., Herbig & Goodrich 1986; 
Herbst et al. 1994; Valcnti, Johns-Krull, & Linsky 2000). 
The UV excess radiation is considered to result from the 
shock on the stellar surface which is caused by a mag- 
netospherical accretion flow from the accretion disk onto 
the star (e.g., Calvet & Gullbring 1998; Ostriker & Shu 
1995). In order to examine the effects of this UV excess 
radiation on the disk structure and the H 2 line emission 
from the disk, we treat two kinds of models for radiation 
from the central star, F„ iSt ar: models with and without 
UV excess. In the model with UV excess, F^ sta r consists 
of three components: black body emission at the star's ef- 
fective temperature, T», an optically thin hydrogenic ther- 
mal bremsstrahlung emission at a higher temperature, and 
Ly a line emission, based on observations towards TW 
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Fig. 2. The vertical temperature profiles of gas (solid lines) and dust (dotted lines) at the disk radii of x = 1, 10 and 
100AU (a-1) (b-1), and the radial temperature profiles of gas at z/x = 0.0,0.1,0.2,0.4, and 0.8 (a-2) (b-2) for the 
models (a) with and (b) without UV excess. In the model (a) the gas temperature at the disk surface can reach around 
1,000 K due to the grain photoelectric heating, while it is almost the same as the dust temperature near the midplanc. 
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Hydrae (see Appendix C). In the model without UV ex- 
cess, the stellar radiation is assumed to be black body 
emission only. 

For the interstellar radiation field, -F^isrf, we adopt 
the Draine (1978) field for the wavelength range of 912A < 
A < 2000A and the field given by van Dishocck & Black 
(1982) for A > 2000A. 



2.3. Dust model 

The physical and chemical structure of protoplanetary 
disks is affected by the dust model in various ways: first, it 
affects the radiation field because dust grains are the dom- 
inant opacity source in protoplanetary disks. Thus, the 
dust temperature profile is influenced by the dust model 
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with UV excess. The grain photoelectric heating dominates the heating process, while the gas-grain collisions and the 
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through the absorption (k^) and extinction (x v ) coeffi- 
cients (Sect. 2.1). Also, it affects the UV radiation field 
(through the extinction coefficient \v and albedo ui v \ Sect. 
2.2), that is, the photodissociation and photoionization 
processes (see Appendix B and Sect. 3). The gas temper- 
ature profile is also affected because grain photoelectric 
heating induced by UV photons is the dominant source 
of gas heating at the disk surface (see Appendix A.l and 
Sect. 2.4) and by grain-gas collisions (see Appendix A. 2). 
Second, the total surface area of the dust grains affect the 
molecular hydrogen abundance through the H2 formation 
rate on grain surfaces (Sect. 3). 

In this paper we assume that the dust particles consist 
of silicate, carboneous grains, and water ice, and have size 
distributions obtained by Weingartner & Draine (2001a; 
hereafter WDOla) which reproduce the observational ex- 
tinction curve of dense clouds (see Appendix D for de- 
tails). Also, the dust and gas are assumed to be well- 
mixed. 



2.4. Resulting temperature and density profiles 

We obtain self-consistent density and temperature dis- 
tributions of a protoplanetary disk by solving the equa- 
tions for vertical hydrostatic equilibrium and local ther- 
mal balance between heating and cooling of gas itera- 
tively together with the local radiative equilibrium equa- 
tions on dust grains (see Sect. 2.1). Figure shows the 
contour plots of the resulting temperature (a-1) and den- 
sity (a-2) profiles in the x-z plane. The contour levels are 
T = 30,100,300, and 1000K for the temperature, and 
p = 10- 20 ,10- 18 ,10- 16 ,10- 14 , and 10~ 12 g cm" 3 for the 
density. In order to see the effects of the UV excess radia- 
tion, we also calculate for comparison the structure of the 
disk which is irradiated by stellar radiation without UV 
excess (Fig.^; see Sect. 2.2). 

Figure |21 shows the corresponding gas (solid lines) and 
dust (dotted lines) temperature profiles in the vertical di- 
rection at the radii of x = 1,10 and 100 AU (a-1) (b- 
1). The horizontal axis is the vertical height divided by 
each radius. The gas temperature profiles as a function of 
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the distance from the central star, R — (x 2 + z 2 ) 1 / 2 , are 
also plotted in (a-2) and (b-2) for z/x = 0.0,0.1,0.2,0.4, 
and 0.8. The figures show that in the model with UV 
excess radiation the gas temperature at the disk surface 
is much higher than the dust temperature and reaches 
around 1,000K even at the radius of 10AU. This is due 
to the grain photoelectric heating induced by the strong 
far ultraviolet (FUV) radiation from the central star. On 
the other hand, near the midplane of the disk and at the 
outer disk in the model without UV excess, where the UV 
radiation field is not strong enough to heat up the gas, the 
gas temperature is close the dust temperature due to the 
collisions between the gas and dust particles (see below). 
We calculate the gas temperature only in the upper re- 
gion of the disk where the difference between the gas and 
dust temperatures are more than 1% (\T — Td\/T > 10 -2 ), 
considering that the gas temperature is almost the same 
as the dust temperature near the midplane. The gas tem- 
perature in ionized region is assumed to be T = 10 4 K 
(see Appendix B.l). Now, the corresponding vertical (1) 
and radial (2) density profiles of the disk in the models 
with (solid lines) and without (dotted lines) UV excess 
radiation are plotted in Fig. The strong UV radiation 
field at the disk surface heats the gas hot enough to make 
the disk expand in the vertical direction. Therefore, the 
higher density at the disk surface in the model with UV 
excess radiation is caused by the stronger UV radiation 
fields. 

The integrated FUV radiation fields (6 eV < hv < 
13 eV) in the disk for the models with and without UV 
excess radiation are shown in Fig. The radial (dashed 
lines), Gfuv,r, and vertical (dotted lines), Gfuv.z, fields 
are calculated by integrating the specific radiation field 
Fv,r (Eq. |B]) and F UjZ (Eq. [J]), respectively, in the 
FUV region (6 eV < hv < 13 eV). The total radiation 
fields (solid lines) are summations of the radial and verti- 
cal fields, Gfuv = Gfuv.r + Gruv.z- The figures show 
that the direct irradiation from the central star domi- 
nates the radiation field in the upper layer, while the 



scattered field is superior in the lower layer (see also van 
Zadelhoff et al. 2003; Bergin et al. 2003). The contribu- 
tion of the interstellar radiation field (dot-dashed lines), 
Gfuv.isrf = / fl.isRF exp[—Tv,z(?<x)]dVi is also plotted 
in this figure, which shows that it is not dominant except 
in the lower layer of the outer region {x ~ 100 AU) in the 
model without UV excess radiation. 

The vertical profiles of the heating and cooling rates 
at the radii of (a) 10AU and (b) 100AU for the model 
with UV excess radiation are plotted in Fig.|SJ We can see 
from the figure that the dominant heating source is grain 
photoelectric heating in the region where we calculate the 
gas temperature (\T - T d \/T > 10~ 2 ). Meanwhile, the 
energy transfer from gas to dust particles via collisions is 
the dominant cooling source at the inner region and near 
the midplane of the disk where the matter is dense enough, 
while radiative cooling via line transition is dominant in 
less dense region. 

3. Abundance and excitation of Ho 

3.1. The model 

The abundance and the level populations of the X X Y>~ 
electronic state of molecular hydrogen in a statistical 
equilibrium state are calculated based on Wagenblast & 
Hartquist (1988) as 



n;(H 2 ) 



Y n m(H 2 ) ( A ml + f3 ml + Y n s c ; 



n(H) J R form , / (8) 



ml. 



where A\ m is the Einstein A-coefficicnt for spontaneous 
emission from level I to level m (Ai m = if the en- 
ergy of level I, Ei, is smaller than E m ) and Cf m is 
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Fig. 7. The level populations of molecular hydrogen at the disk radii of x = 1 AU (filled diamonds), 10AU (open 
squares) and 100AU (asterisks) for the models (a) with and (b) without UV excess radiation. The LTE distributions 
are plotted in dashed lines. If the central star has a UV excess, the level populations are in LTE at the inner disk, and 
the populations of the upper levels become large due to the high temperature. 
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the collisional transition rate with collision partner s. 
The collisional rate coefficients are taken from Smith 
et al. (1982) and Black & Dalgarno (1977) for colli- 
sions with H + , Tine et al. (1997; see also references 
therein) for the collisions with H 2 , and Lepp et al. (1995; 
see also Tine et al. 1997) for the collisions with H 
( [http://www.physics.un lv.edu/astrophysics / ) . 

The symbol j3i m in Eq. JSJ) represents the effective rate 
for transition I — > m via ultraviolet pumping followed 
by radiative cascades and i?diss,z is the photodissociation 
rates of hydrogen molecules in level The UV radiation 
fields are calculated as discussed in Sect. 2.2. An approxi- 
mate function given by Federman et al. (1979) is used for 
the H 2 self-shielding of the UV radiation fields. A purely 
thermal Doppler widths of Vd = (2fcT/m Al ) 5 is adopted 
in this function. Molecular hydrogen excitation due to X- 
ray induced electron impact should be taken into account 
in future (e.g., Bergin et al. 2004). 

The symbol R{ or m,i represents the effective formation 
rate of H2 in level I on grain surfaces. The molecular hy- 
drogen which leaves the grain surfaces is assumed to orig- 
inate in the levels (v — 7, J = 1) and (v = 7, J = 0) 
with the ratio of 9:1 and cascade into the level I, following 
Duley & Williams (1986) (cf. Takahashi & Uehara 2001 
for new formation pumping model). The total formation 
rate of J2i flform.i = 7.5 x 10- 18 T°- 5 e H2 n tot n(ff) cm' 3 s" 1 
is adopted here. The symbol en 2 represents the recombi- 
nation efficiency of atomic hydrogen on dust grains which 
is estimated by Cazux & Tielens (2002, 2004) based on 
laboratory experiments (e.g., Pirronello et al. 1999; Zecho 
et al. 2002). 

The endothermic reaction O + H2 — > OH + H, which 
dominates the destruction of molecular hydrogen in high 
temperature regions (e.g., Storzer & Hollenbach 1998), is 
also taken into account. The reaction rate coefficient of 
fc 0+H2 = 3.14 x 10- 13 (T/300X)exp(-3150X/T) cm 3 s" 1 
is adopted from the UMIST RATE99 database (Le Teuff 
et al. 2000). 

The total number density of hydrogen nuclei at each 
point in the disk is required to satisfy the condition, 

n tot = n(H) + 2n(H 2 ) + n(H+) + 2n(H+) + 3n(H+), (9) 

where n(s) represents the number density of species s. A 
simple chemical equilibrium scheme given in Wagenblast 
& Hartquist (1988) is used to obtain the number densities 
of the related species. 

3.2. Results 

The resulting molecular (solid lines) and atomic (dashed 
lines) hydrogen abundances in the vertical direction at 
the radii of x = 1, 10 and 100 AU are plotted in Fig. 
Molecular hydrogen is photodissociated at the very surface 
of the disk, while it is destroyed by atomic oxygen where 
the gas temperature is ~ 1, 000 K. 

Figure[7]shows the resulting level populations of molec- 
ular hydrogen. The filled diamonds, open squares, and as- 
terisks show the column densities of molecular hydrogen in 



each ro- vibrational level as a function of the level energy at 
each radius of 1, 10 and 100 AU, respectively. The column 
densities are calculated by integrating the number density 
of molecular hydrogen in each level along the vertical di- 
rection at each radius of the disk. The level populations 
in local thermodynamic equilibrium (LTE) at each radius 
are shown in the dashed lines. Also, we present the level 
populations in the model without UV excess in Fig. [TJs 
for comparison. Figures labeled (1), (2), and (3) show the 
level populations at the radius of 1AU, 10AU, and 100AU, 
respectively, and they are plotted together in figure (0) for 
comparison. We can see from the figures that if the cen- 
tral star has a UV excess, the gas becomes hot enough, as 
we have seen in the previous section, that the collisional 
processes are very efficient and the level populations are 
in LTE except for those of the very high levels. And as 
a result of the high temperature, the populations of the 
upper levels become very large. On the other hand, in the 
outer disk in the model without UV excess where the gas 
is cold enough, the level populations are not in LTE but 
are affected by the UV pumping process. Also we can see 
that in this model, the populations of the upper levels are 
not as large as those in the model with UV excess. 

4. Molecular hydrogen emission 

Making use of the level populations we obtained in Sect. 
3, we calculate the molecular hydrogen emission from the 
disk. The intensity of each line integrated over the fre- 
quency, I u i, is obtained by solving the radiative transfer 
equation in the vertical direction of the disk, assuming 
that an observer faces the disk, 

-T^ = ~Xul{Iul - S u l), (10) 

az 

where the subscripts ul means the transition from the up- 
per to the lower levels. The source function, S u i, and the 
total extinction coefficient, \ui, are given as 

S ut = —n u A ul <i> u i^ (11) 
Xui 4tt 

and 

Xui = PXv ul + {niB lu - n u B u i)^ u i^^-, (12) 

where the symbols A u i and B u i are the Einstein coeffi- 
cients for the transition u — > I, and n u and ni are the 
number densities of the upper and lower levels, respec- 
tively. The energy difference between the levels u and I 
corresponds to hv u i- The symbol is the line profile 
function, p the gas density, and Xv u i the dust opacity at 
the frequency v u i- 

The observable line flux is obtained by integrating Eq. 
(|10fl in the vertical direction and summing up the integrals 
in the radial direction of the disk as, 

Ful = ——2 I 2nxdx / rj u i(x,z)dz, (13) 
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Fig. 8. The near-infrared (1/im < A < 4/im) line spectra of molecular hydrogen from the disk in the model (a) with 
and (b) without UV excess radiation of the central star. 
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Fig. 9. Same as Fig. [SJ but for the mid-infrared (5^m < A < 30/im) wavelength band. 



where fj u i(x 7 z) is the emissivity of the transition line at 
(x, z) times the effect of absorption in the upper disk layer, 

rj u i(x,z) = n u (x,z)A u i^^exp(-T u i(x, z)), (14) 

and t u i(x, z) is the optical depth from z to the disk surface 
ZoQ at the frequency v u i, 

t u i(x,z) = J Xul{x, z')dz' '. (15) 

Here, we use the distance to an object of d = 56 pc for 
calculating the intensity in order to compare it with the 
observations towards TW Hya. 

4.1. The infrared spectra 

The resulting line spectra in the near- and mid-infrared 
wavelength bands are plotted in Figs. IHJi and|5^, respec- 
tively. The line spectra for the model without UV excess 
are presented in Figs.|Sb and|5jD for comparison. We can 
see from the figures that if the central star has a UV ex- 
cess, the line spectra become stronger by about a few or- 
ders of magnitude because the strong UV irradiation by 
the central star heats the gas hot enough and makes the 



populations of the upper levels large, as we have seen in 
Sect. 3.2. 

Some of the strongest line fluxes (Fli ne ) we calculate 
for the models with (Model A) and without (Model B) UV 
excess radiation are listed in Tabled The intensity ratios 
of the 2.Yl[im v = 1 -> 5(1) line to the 2.25/xm v = 
2^1 5(1) line are ~ 80 and - 18 for Model A and B, 
respectively. The corresponding excitation temperatures 
are ~ 900 K and ~ 1300 K, which are comparable to the 
kinetic temperature (cf. Fig. |3J) as the level populations 
are in LTE (see Fig. 0. The excitation temperature is 
higher in the model without UV excess than in the model 
with UV excess because the emission comes mainly from 
the inner disk in the former case while it comes mostly 
from the outer cooler disk in the latter case as we will 
see below. The line ratio of the 2.12/^to v = 1 — > 5(1) 
line to the 2.41/ito v = 1 — > Q(l), two of the strongest 
lines in Fig. [SI is larger for the model with UV excess for 
the same reason. On the other hand, the pure rotational 
transition lines in Fig. ED are emitted from all over the 
disk for both models since they trace the gas in cooler 
regions than the vib-rotational transition lines in Fig. [S] 
(see below). Therefore, the excitation temperature derived 
from the pure rotational transition lines for the model with 
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Table 1. The observed and calculated infrared line flux of molecular hydrogen (-Fime), the calculated continuum 
radiation flux (F cont ) [erg/s/cm 2 ], and the line-to-continuum flux ratio (Fy mc / F cont ). 



X 




Line 


Model A 


^linc 

Model B 


Observation 


Tcont . 




(/im) 






(xlO" 15 ) 


(xlO" 18 ) 


(xlO" 15 ) 


(xlO" 9 ) 


(xlO" 6 ) 


1.16 


2 


-0 S(l) 


0.03 


0.10 




5.19 


0.01 


1.19 


2 


- 5(0) 


0.01 


0.05 




5.09 


0.00 


1.24 


2 


-0 Q(l) 


0.03 


0.44 




4.91 


0.01 


1.24 


2 


-0 Q(2) 


0.01 


0.05 




4.90 


0.00 


1.25 


2 


- Q(3) 


0.02 


0.06 




4.87 


0.00 


1.83 


1 


- 5(5) 


0.40 


0.13 




2.87 


0.14 


1.89 


1 


- S(4) 


0.31 


0.11 




2.73 


0.11 


1.96 


1 


- S(3) 


1.81 


0.70 




2.57 


0.71 


2.03 


1 


- S(2) 


0.95 


0.47 




2.40 


0.39 


2.12 


1 


-0 S(l) 


3.43 


2.44 


1.0 a 


2.22 


1.55 


2.22 


1 


- 5(0) 


0.91 


0.89 




2.03 


0.45 


2.25 


2 


- 1 S(l) 


0.04 


0.13 




1.99 


0.02 


2.35 


2 


- 1 S(0) 


0.01 


0.07 




1.81 


0.01 


2.41 


1 


-0 Q(l) 


3.56 


7.09 




1.74 


2.05 


2.41 


1 


-0 Q(2) 


1.00 


1.00 




1.73 


0.58 


2.42 


1 


- Q(3) 


2.41 


1.78 




1.71 


1.41 


2.44 


1 


- Q(4) 


0.53 


0.27 




1.69 


0.31 


2.45 


1 


-0 Q(5) 


0.88 


0.34 




1.67 


0.53 


2.47 


1 


- Q(6) 


0.14 


0.05 




1.64 


0.09 


2.50 


1 


- Q(7) 


0.18 


0.06 




1.61 


0.11 


2.55 


2 


- 1 Q(l) 


0.05 


0.69 




1.55 


0.03 


2.56 


2 


- 1 Q(2) 


0.01 


0.08 




1.54 


0.01 


2.57 


2 


- 1 Q(3) 


0.03 


0.10 




1.52 


0.02 


6.11 





- 3(6) 


1.22 


0.62 




0.25 


4.98 


6.91 





-0 5(5) 


6.25 


4.14 




0.19 


32.06 


8.02 





- 5(4) 


2.90 


2.67 




0.18 


16.10 


9.66 





- 5(3) 


9.26 


12.62 




0.24 


38.37 


12.27 





- 5(2) 


2.17 


5.33 


< 30 c 


0.15 


14.28 


17.02 





- 5(1) 


2.38 


14.33 


28-8l\< 39 c 


0.13 


18.62 


28.20 





- 5(0) 


0.08 


3.08 


25 - 57 b 


0.09 


0.94 



a Observation by Bary et al. (2003). 

b Observations by Thi et al. (2001b). 

c Observations by Richter et al. (2002). 

d Fi in e in Model A. 



UV excess is higher (that is, the weaker line flux for the 
lower level transitions) than that derived from the lines for 
the model without UV excess because the UV irradiation 
from the central star heats the gas in the disk more in the 
former model. 

The observational intensity of the ro-vibrational line, 
v = 1 -> 5(1), towards TW Hya by Bary et al. (2003) 
is also presented in Tabled an d is comparable to our cal- 
culated flux in the model with UV excess. In addition, 
the observations of mid-infrared line spectra of pure ro- 
tational transitions, 5(0) and 5(1), by ISO (Thi et al. 
2001b) and the upper limit of the line fluxes of 5(1) and 
5(2) under the assumption of FWHM of 30 km s" 1 ob- 
served by TEXES at the NASA IRTF (Richter et al. 2002) 
are presented in Tabled These are observed towards clas- 
sical T Tauri stars in the Taurus-Auriga cloud complex 
(d = 140pc). We can see from the table that the calculated 



line fluxes are weaker than the results of the ISO observa- 
tions, and consistent with the ground-based observations 
by TEXES. It may be that the spatial and wavelength 
resolutions of the ISO observations were not high enough 
so that the observed line emission comes not from the cir- 
cumstellar disk but from shocked regions in the associated 
outflows. The last two columns of Table H represent the 
flux of the calculated dust continuum emission from the 
disk plus the radiation from the central star (-Fcont.) and 
the flux ratio of the molecular hydrogen line in Model A 
to the continuum emission (i^i no /-Pcont.)- The dust con- 
tinuum emission is calculated by solving the radiative 
transfer equation (Eq. ^U|), adopting S u i = B Vul (T ( i) and 
Xui = PXu^i for the source term and the opacity, respec- 
tively. The central star is assumed to emit blackbody ra- 
diation with an effective temperature of T» = 4000K and 
a radius of R* = 2Rq. 
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Fig. 10. The radial flux distributions of the 2.12/im v = 
1 — 5(1) line and the 28/im 5(0) line for the models with 
(solid line) and without (dashed line) UV excess radiation 
(top). The maximum gas temperature at each radius is 
also plotted at the bottom. If the central star has UV 
excess radiation, the 2.12/im v = 1 — 5(1) line emis- 
sion comes mainly from the radius of 20 AU, while the 
28/im 5(0) line emission comes from all over the disk. 



Figure 1101 represents the radial flux distributions, 
dF u i/dx = {2-kx/ Ai:d 2 ) J 00 fj u idz, of the 2.12/im v = 
1 — ► 5(1) line and the 28/im 5(0) line for the mod- 
els with (solid line) and without (dashed line) UV ex- 
cess radiation. The maximum gas temperature at each 
radius, T maK (x) — max{T(a;, z)|0 < z < Zoo}, is also 
plotted at the bottom of the figure. We can see from 
the figure that if the central star has a UV excess, the 
2. 12 fim v = 1 — > 5(1) line emission comes mainly from 
around the radius of 20 AU, which is consistent with that 
derived from the observed peak-to-peak velocity separa- 
tion of the H2 emission line due to the rotation veloc- 
ity of the disk (Bary et al. 2003) . The line flux distribu- 
tion increases with increasing radius inside 20 AU since 
the line emissivity integrated over the vertical direction, 
f_°° Vuidz, does not change with radius and dF u i I dx is al- 
most proportional to the radius. This is because the tran- 
sition line arises from hot gas with T > 1000K and the 
maximum gas temperature of the disk satisfies this con- 
dition. The maximum temperature does not change very 
much at x < 20AU due the weak dependence of the photo- 
electric heating rate on the UV photon flux in this range 
(see e.g., Weingartner & Draine 2001b). Meanwhile, the 
line flux distribution decreases suddenly beyond 20 AU be- 
cause the gas temperature falls rapidly below 1000 K due 
to the stronger dependence of the photoelectric heating 
rate on the UV flux in this range. For the model without 
UV excess, the line emission mainly comes from the very 
inner region of the disk, x < 0.3AU, because the UV irra- 
diation by the central star is not strong enough to heat the 
gas at the outer disk to high enough temperatures as we 
have seen in Sect. 2.4. The line flux distribution increases 
with increasing radius at the outer disk, x > 3AU, since 
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Fig. 11. The vertical profiles of the emissivity of the 
2.12/xm v = 1 - 5(1) line and the 28/im 5(0) line at 
the radius of 10 AU for the models with (solid line) and 
without (dashed line) UV excess radiation (top). The ver- 
tical profiles of the gas temperature and density are also 
plotted at the bottom. For the model with UV excess radi- 
ation, the 2.12/im v — 1 — 5(1) line emission comes from 
hot surface layer with T gas > 1000K, while the 28/im 5(0) 
line emission comes from cooler region near the midplane 
with T gas > 50K. 



the populations of the ro-vibrational levels are controlled 
not by the collisional excitation but by the UV pumping 
process (see Fig.^t 1 ), and the radiation energy density in- 
tegrated over the vertical direction is almost constant to 
the radius in this region. 

The 28/im 5(0) line emission comes from a cooler re- 
gion, compared with the 2.12/im v — 1 — > 5(1) line, with 
the gas temperature of T > 100 K. If the central star has 
UV excess radiation, the radial flux distribution of the line 
is almost constant all over the disk, that is, the line emis- 
sivity integrated over the vertical direction decreases with 
increasing radius. This is because the gas is heated up to 
T > 100 K even in the dense region near the midplane in 
the inner disk, while it is heated up only in the less dense 
region of the disk surface in the outer disk (see Fig. [SJi) . 
Meanwhile, for the model without UV excess, the gas tem- 
perature does not reach 100 K at x > 20 AU (see Fig. 
which makes the line flux very weak at the outer disk. The 
flux distributions for the models with and without UV ex- 
cess are identical at the inner disk, x < 3AU, because the 
gas temperature reaches 100 K even in the region where 
it is controlled by the dust temperature through gas-grain 
collisions, and the dust temperatures are identical in both 
models (see Fig. - 

Furthermore, the vertical profiles of the emissivity, 
rjul (Eq. H3J), of the 2.12/im v = 1 -> 5(1) line and 
the 28/im 5(0) line at 10 AU for the models with (solid 
line) and without (dashed line) UV excess radiation are 
plotted in Fig. [TT] We can see from the figure that in 
the surface layer of the disk the emissivity of both lines 
for both models basically increases with decreasing ver- 
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Table 2. The observed and calculated ultraviolet line flux of molecular hydrogen [erg/s/cm 2 ] 



A 




Line 


Model A 


Model B Model C 


Obs. a 


A 




Line 


Model A 


Model B Model C 


Obs. a 


(A) 






(xKT 15 ) 


(xl0~ 21 


) (xl0~ 16 ) 


(xKT 15 ) 


(A) 






(xlO" 15 ) 


(xl0~ 21 


) (xl0~ 16 ) 


(xlO- 15 ) 


1161.7 





- 1 P(0) 


0.16 


4.42 


0.08 




1148.7 


1 


- 1 R(3) 


0.65 


0.04 


0.18 


4.6 


1166.3 





- 1 -P(2) 


0.25 


8.74 


0.14 




1161.9 


1 


- 1 -P(5) 


1.25 


0.05 


0.30 


10.9 


1217.3 





- 2 R(0) 


2.61 


13.22 


0.37 




1202.5 


1 


- 2 i?(3) 


7.04 


0.08 


0.79 


11.3 


1222.0 





- 2 P(2) 


4.40 


26.14 


0.71 




1216.2 


1 


- 2 P(5) 


13.11 


0.09 


1.15 




1274.6 





- 3 R(0) 


7.75 


23.68 


0.70 


27.4 


1257.9 


1 


- 3 i?(3) 


17.68 


0.06 


0.89 


18.1 


1279.6 





- 3 P(2) 


14.86 


46.82 


1.37 


39.2 


1272.0 


1 


-3 P(5) 


23.49 


0.07 


1.12 


20.5 


1333.6 





- 4 R(0) 


9.64 


28.21 


0.84 


42.8 


1314.8 


1 


- 4 R(3) 


3.81 


0.01 


0.17 


12.2 


1338.7 





-4 P(2) 


18.83 


55.75 


1.65 


73.1 


1329.3 


1 


-4 P(5) 


4.61 


0.01 


0.21 


7.5 


1393.9 





- 5 R(0) 


7.94 


23.21 


0.69 


35.3 


1372.7 


1 


- 5 R(3) 


2.67 


0.01 


0.12 


3.2 


1399.1 





- 5 P( 2 ) 


15.51 


45.88 


1.36 


73.8 


1387.5 


1 


- 5 P(5) 


3.23 


0.01 
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Observations by Herczeg et al. (2002). 



tical height due to the increasing gas density. For the 
model with UV excess radiation, the emissivity of the 
2.12yum v = 1 — » 5(1) line decreases suddenly where the 
gas temperature falls below ~ 1000 K. The correspond- 
ing disk height is z ~ 3.5AU^ 17 H, where H = c s /£l is 
the disk scale height. The optical depth from the disk sur- 
face to z ~ 3.5AU is small enough, t u i ~ 1.5 x 10~ 3 , for 
the line. Meanwhile, in the model without UV excess, it 
is independent of the gas temperature because the level 
populations are controlled by the UV pumping process 
(see Sect. 3.2), and it decreases near the midplane where 
the UV irradiation cannot penetrate. The emissivity of 



the 28/zto 5(0) line has a peak where the gas temperature 
drops below 100K in both models. The corresponding disk 
height is z <~ 1.9AU^ 9.5H, where the line is also optically 
thin (t u i ~ 6.6 x 10~ 3 ) for the model with UV excess. The 
emissivity in the far cooler region close to the midplane 
is still relatively large due to the high gas density, but 
it decreases distinctly at T < 50 K. Therefore, although 
the gas in the disk is concentrated near the midplane, 
these lines come only from the upper region of the disk. 
Also, the figure indicates that the line flux traces the gas 
density and temperature where the gas temperature falls 
below the critical value. 
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Fig. 13. The ultraviolet line spectra of molecular hydro- 
gen for the model without Lya line emission from the 
central star. 
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Fig. 14. The level populations in the excited electronic 
state at (X, Z) = (10AU, 6AU) for the models with UV ex- 
cess (asterisks) and without Lya emission (open squares). 
Some specific levels are highly populated due to irradia- 
tion by the strong Lya line emission from the central star. 



4.2. The ultraviolet spectra 

The resulting line spectra in the ultraviolet wavelength 
band are plotted in Figs. 112b and 1 12b for the models with 
and without UV excess radiation, respectively. If the cen- 
tral star has a UV excess, some of those lines, whose upper 



levels are in excited electronic states and populated via the 
UV pumping process, become stronger by more than six 
orders of magnitude than those in the model without UV 
excess. 

Some of the strongest line fluxes we calculate in the 
model with UV excess (Model A) are listed in Tableland 
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1 — > 3 P(5) line for the models with (solid line) and with- 
out (dashed line) UV excess radiation, and without Lya 
emission (dotted line). The radiation flux decreases with 
increasing radius for all models because the UV irradia- 
tion flux, which pumps hydrogen molecules in the ground 
electronic state, decreases. 
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Fig. 16. The vertical profiles of the emissivity of the 
1272.0A t!=l->3 P(5) line at the radius of 10 AU 
for the models with (solid line) and without (dashed line) 
UV excess radiation, and without Lya emission (dotted 
line) (top). The energy density profiles are similar to the 
number density profiles of H2 in the pre-excited level in 
the ground electronic state (bottom). 



compared with observations by the HST and the FUSE 
(Herczeg et al. 2002). The calculated strong line fluxes re- 
sult from the excitation due to the strong Lya line emis- 
sion from the central star and are consistent with the ob- 
servations. The line fluxes in the model without UV excess 
(Model B) are also presented in the table for comparison. 
In addition, in order to see the effect of pumping by Lya 
emission, we calculate the level populations and UV line 
fluxes by neglecting Lya emission (Model C). The gas den- 
sity and temperature profiles for the model with UV excess 
are used here. The resulting line emission is presented in 
Fig-ED an d Tabled and is much weaker than the observed 
values. 

The level populations in the excited electronic state at 
the radius of 10 AU and the vertical height of 6 AU for 
Model A (asterisks) and C (open squares) are plotted in 
Fig-El Figures labeled (1) and (2) show the level popula- 
tions for Model A and C, respectively, and they are plotted 
together in figure (0) for comparison. The dotted lines in 
figure (0) connect the level populations of Model A and 
C in the same energy levels. The figures clearly show that 
specific level populations are extremely high for Model A 
(with UV excess radiation) due to the irradiation of the 
Lya emission from the central star. The transitions from 
these levels result in the strong line emission in Fig. 112k . 
Meanwhile, the level populations for Model C (without 
Lya emission) are low and distributed smoothly since they 
are excited only by continuum radiation. This is why the 
line fluxes are weak and many lines have similar strength 

in Fig. na 

In the case of Model B (without UV excess) the line 
spectra are much weaker (Fig. 112b ) due to the weak UV 
irradiation from the central star and the low level popula- 
tions of pre-excited hydrogen molecules in the ground elec- 



tronic state (Fig. 03, see also below). Also, the line spec- 
tra are relatively scattered, although they result from the 
pumping by continuum irradiation. This is because molec- 
ular hydrogen in highly excited electronic states originates 
mainly from highly populated low energy levels (v = 0) 
in the ground electronic state, whose populations decrease 
suddenly with increasing level energies due to the low tem- 
perature for the model without UV excess (Fig. 0>)- Thus, 
the level populations in the excited electronic state also 
have steeper distribution, which leads to the sparse line 
spectra. 

Figure shows the radial flux distribution of the 
1272.0A v = 1 -> 3 P(5) line, which is one of the strongest 
lines in Model A, for the models with (solid line) and with- 
out (dashed line) UV excess radiation, and without Lya 
emission (dotted line) . The flux decreases with increasing 
radius for all models because the line intensity is propor- 
tional to the population of the upper transition level in the 
excited electronic state, namely, the UV irradiation flux 
from the central star which is proportional to the inverse 
squares of the radius at the disk surface. 

The vertical profiles of the line emissivity at 10AU for 
these models are also plotted in Fig. ED (top). At the 
disk surface, the profiles of f) u i are similar to those of 
the number densities of pre-excited hydrogen molecules in 
the ground electronic state, n(v, J), plotted in the bottom 
of Fig. 1161 This is because the line emissivity is propor- 
tional to the number density of H2 in the upper tran- 
sition level of the line, n(v* = I, J* = 4) (the sym- 
bol * means that it is in the excited electronic state), 
which is populated via the UV pumping process and pro- 
portional to both n(v, J) and the UV irradiation flux, 
K(vJ^v* j*),r + Fv(vj->vj*),z (Eqs. |B] and 0). In the 
models without UV excess and without Lya emission, 
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the level (v* = 1, J* = 4) is populated mainly from 
the highly populated low energy level, (v = 0, J = 3). 
Meanwhile, in the model with UV excess, the transition 
from the level (v = 2, J = 5) dominates in populating the 
level (v* = 1, J* =4) because the excitation wavelength, 
A = hc/E(v = 2, J = 5 -»• v* = 1, J* = 4) = 1216.2A, is 
close to the central wavelength of the Lya line emission. 
So, the number density of n(v = 0, J = 3) is plotted for 
the models without UV excess and without Lya emission, 
while that of n(v = 2, J = 5) is plotted for the model with 
UV excess in the bottom of Fig. QUI The line emissivity 
decreases with decreasing vertical height for all models 
especially near the midplane because the irradiating UV 
photons are self-absorbed (used for the pumping) and can 
not penetrate deep in the disk. Thus, the molecular hy- 
drogen emission from near the midplane is not observable 
in the ultraviolet region. 

5. Summary 

We have, for the first time, investigated molecular hydro- 
gen emission from protoplanetary disks, taking into ac- 
count the global physical structure of the disk and the 
detailed analysis of level populations of molecular hydro- 
gen. 

First, we obtained in a self-consistent manner the den- 
sity and temperature profiles of gas and dust, taking into 
account the irradiation by a central T Tauri star. As a re- 
sult, we found that if the central star has UV excess emis- 
sion over the photospheric blackbody radiation, the gas 
temperature at the disk surface can reach about 1,000K 
even at the radius of 10 AU due to grain photoelectric 
heating induced by the UV photons from the central star. 

Next, making use of the physical structure of the disk, 
we calculated the level populations of molecular hydrogen 
in its ground electronic state on the assumption of statis- 
tical equilibrium. The resulting level populations tend to 
be in LTE, although the populations in the upper levels 
become very large due to the high gas temperature if the 
central star has UV excess radiation. 

Furthermore, using these level populations, we com- 
puted molecular hydrogen emission in the near- and mid- 
infrared and ultraviolet wavelength bands by solving the 
radiative transfer equation. Consequently, we found strong 
line emission spectra if the central star has UV excess ra- 
diation. The infrared line spectra are strong because the 
populations of the upper levels become large due to the 
high gas temperature, while the ultraviolet line emission 
is intense because the strong UV irradiation, including 
that due to Lya, by the central star pumps the hydrogen 
molecule from the ground to the excited electronic states. 
We compared the results of our calculation with the obser- 
vations to find that the calculated near-infrared intensity 
of the u = l— ►OiS'(l) line is in good agreement with that 
observed towards TW Hya, as are the ultraviolet lines. 
Our predictions for the pure rotational line intensities are 
much below those observed by Thi et al. (2001b) but this 
may be related to the large ISO beam used. 
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Appendix A: Thermal processes 

The following three processes which are dominant under 
conditions in protoplanetary disks are considered as heat- 
ing and cooling processes in this paper. 

A.l. Grain photoelectric heating 

The photoelectric emission from dust grains induced by 
far ultraviolet (FUV) photons from the central star dom- 
inates the heating at the surface layer of protoplanetary 
disks (see Sect. 2.4), especially if the star has UV excess ra- 
diation (Sect. 2.2 and Appendix C). In this paper we adopt 
the photoelectric heating rate calculated by Weingartner 
& Draine (2001b), 

r pc = 1.0 x 10~ 26 eG F uvntotergs cm" 3 s _1 , (A.l) 



1 + G 2 (G FU v VT/n e f* [1 + G 3 (G FUV VT/n e f°} ' 

where Gfuv represents the FUV fields measured in unit 
of the equivalent average interstellar radiation flux of 
1.6 x 10 -3 ergs cm" 2 s" 1 (Habing 1968), and T, n tot , 
and n e represent the gas temperature, the number den- 
sities of hydrogen nuclei and electrons, respectively. The 
parameter set of Go - C§ is taken from Weingartner & 
Draine (2001b), where we choose the dust size distribution 
model for dense clouds with the ratio of visual extinction 
to reddening of R v = A(V)/E(B -V) = 5.5, the total C 
abundance per H nucleus in the log-normal population of 
6c = 3.0 x 10~ 5 , and the grain volumes assumed to be the 
same as those of the diffuse interstellar medium (Case B 
of WDOla; see also Appendix D). 

A. 2. Gas-grain collisions 

The energy exchange between gas and dust particles 
through collisions is one of the most important thermal 
processes in protoplanetary disks in which the density is 
high enough. We use the following cooling rate (heating 
rate if Td > T) in this paper, 

A gr = 3.5 x 10- 34 n 2 ot T a5 (T - T d )ergs cnrV 1 (A.3) 

(Burke & Hollenbach 1983; Tielens & Hollenbach 1985). 
The dust temperature Td is obtained on the assumption 
of the local radiative equilibrium in Sect. 2.1. 
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A. 3. Radiative cooling 

Radiative transitions among the fine-structure levels of 01 
(63^m) and CII (158/xm), and the rotational line transi- 
tions of CO contribute to cooling gas in the surface layer 
of the disks where the gas density is low enough. The cool- 
ing rate due to the transition from upper level i to lower 
level j of species s is calculated as 

A a (i/y) = hvijP eS c(Tij) 



x //,{.!,, + />',,/';(',, )} - TijBjiPivij)}, 



(A.4) 



where n;, A 



B-ij, and hv^ are the population density 



of level i, the Einstein probabilities for spontaneous and 
stimulated emission, and the energy difference between 
levels i and j, respectively (e.g., de Jong et al. 1980). 
As the photon escape probability, /? e sc(fij), we use the 
approximate function given by de Jong et al. (1980). 
Assuming that the nearest boundary from a point (r, z) 
is the disk surface (r, Zoo), we approximately calculate the 
optical depth averaged over the line, , as 
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(A.5) 



(A.6) 



where <?i is the statistical weight factors of the level i, 
5vd{z) = {2kT(z) j fim} - 5 is the Doppler line width, and 
Nn(z) = J °° n-f[(z')dz' is the column density of hydrogen 
nuclei from z to z^. The background radiation P(vij) is 
calculated by assuming that it is contributed by the radi- 
ation from the central star and the infrared dust emission, 



P(vij) = (l/8)(i?/i?0 2 5(t/y,T*)exp[ 
+ (l/2)B(z/ y -,T d ){l-exp[- 



-Td,R{Vij)\ 
'*)]}. 



Td.z{ V K 



(A.7) 



where T* and i?* are the temperature and the radius of 
the central star, respectively, and Td t R(vij) — J R n Uij pdR 
and Td, z {vij) = f z °° n Vij pdz are the optical depth in the 
radial and the vertical directions, respectively. The dust 
temperature Td is obtained on the assumption of local 
radiative equilibrium, as discussed in Sect. 2.1. 

The level populations of CII and 01 are obtained by 
solving the equations of statistical equilibrium, 



where 



(A. 



Appendix B: Chemistry 

B.l. Hydrogen photoionization 

It is known that extreme UV photons with energy of 
hv > 13.6 eV are mainly consumed to photoionize atomic 
hydrogen. Thus, we simply consider that most hydrogen 
atoms are ionized if 



a H n tot dr < 

R* J ^EUV 



F, 



v,R 

hv 



dv. 



(B.l) 



where qh = 2.58 x 10~ 13 cm 3 s _1 is the recombina- 
tion coefficient for atomic hydrogen at 10 4 K, F Vt u is 
the radiation flux from the central star in Eq. ©, and 
hvEuv = 13-6 eV (e.g., Osterbrock 1989; Hollenbach et 
al. 1994). The temperature of this ionized region is simply 
assumed to be T — 10 4 K, since we are interested in the 
structure of non-ionized region and the temperature of this 
ionized surface layer does not affect the global structure 
of the disks very much. 

B.2. Carbon and oxygen chemistry 

In order to calculate the radiative cooling by C + , O and 
CO, we treat a very simple chemistry for carbon and oxy- 
gen, and calculate the abundances of these species. The 
total fractional abundances of carbon and oxygen with re- 
spect to hydrogen nuclei are assumed to be a;(Ctot) — 
x(C+) + x(C) -I- x(CO) = 7.86 x 10~ 5 and x(O tot ) = 
x(O) + x(CO) = 1.8 x 10~ 4 , respectively. The number 
density of carbon monoxide is calculated using the follow- 
ing two chemical equilibrium: 

fc n(C+)n(H 2 ) = fcin(OT B )n(0)+GWrcH.n(CH 3 )(B.2) 
and 



fcin(CrL>(0) = G F uvr C on(CO), 
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(A.9) 



i '.7 

HJ — j3ij J 

The symbol represents 
rate. The collisional de-excitation rate coefficients are 
taken from Tielens & Hollenbach (1985), and the 



i > 3, 
i < j. 

the collisional transition 



following Nelson & Langer (1997) (cf. Langer 1976), where 
ko = 5 x 10~ 16 cm 3 s _1 is the rate coefficient of the 
radiative association reaction C + +H2 — >CHj + hv and 
ki = 5 x 10~ 10 cm 3 s _1 is the rate coefficient of the reac- 
tion between the hydrocarbon radicals and atomic oxygen 
to form CO. The symbol Gfuv represents the UV radia- 
collisional excitation rates are calculated as Cy = tion described in Sect. 2.2, and Tch,.. = 5 x 10~ 10 s _1 and 
Cji(gj/gi)exp(huij/kT)(i < j). The level population of 
CO is assumed to be in LTE, which will overestimate the 



r co = io 



-10 



CO cooling rate when nn < 10 cm (cf. Kamp k, van 
Zaldelhoff 2001). The abundances of OI, CII, and CO are 
calculated as is described in Appendix B.2. 



s _1 are the photodissociation rates of the 
hydrocarbon radicals and the carbon monoxide, respec- 
tively. More detailed chemistry should be solved in future 
(e.g., Aikawa et al. 2002; Markwick et al. 2002; Kamp & 
Dullemond 2004). 
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Appendix C: Central stellar radiation 

Here we model the radiation from the T Tauri star 
by means of black body emission at an effective stel- 
lar temperature plus optically thin hydrogenic thermal 
brcmsstrahlung emission at a higher temperature (e.g., 
Lago et al. 1984; Costa et al. 2000) and Ly a line emis- 
sion (e.g., Hergzeg et al. 2002). The specific luminosity at 
a wavelength A of the photosphcric black body radiation 
is calculated as 



10 



10 



^A,bb 



4ir 2 RiB x {T* 



(CI) 



while the luminosity of the coronal hydrogenic thermal 
bremsstrahlung emission is computed as 



(C.2) 



where V is the volume of the bremsstrahlung emission 
region and r?A,br is the emissivity, 



r, x , hl = 2.0 x 10- 36 nX 1/2 A- 2 



x exp(— hc/XkT hr ) ergs 1 cm 3 A 



(C.3) 



(e.g., Rybicki & Lightman 1979). Here n e , Tb r , and c 
are the electron number density, the temperature of the 
bremsstrahlung emission, and the light speed, respectively. 
In addition, we take into account the luminosity of the Ly 
a line emission at the stellar surface, 



L\.Lya — Lx ,Lya exp{ — (A - \ ) 2 / <T 2 } , 



(C.4) 



where we simply assume that the line has a Gaussian pro- 
file with the central wavelength of Ao, the peak luminosity 
of I<Ao,Lya, and the line width of a. 

Now, we choose the physical parameters in Eqs. I|C.2|) - 
l|C.4|l by comparing the calculated radiation flux to obser- 
vations towards TW Hya. The effective stellar tempera- 
ture and radius of TW Hya of T* = 4000K and i?„ = li? , 
and the distance to TW Hya of d = 56 pc are adopted in 
the calculation. Figure IC"T1 shows the calculated radiation 
flux density and the observational data. The dashed, dot- 
ted, dot-dashed, and solid lines represent the black body, 
the bremsstrahlung, the Ly a line, and the total radiation, 
respectively. The filled diamonds are the observations. The 
observational data at the I, R, V, B, and U bands are taken 
from Herbst et al. (1994). The median value of 11.02 mag 
is adopted for the V magnitude. The data of 2.0 xl0~ 14 
and 1.0 xl0~ 14 erg cm~ 2 s~ 4 A -1 at the wavelengths of 
1300A and 3000A, respectively, are measured from the ob- 
servation by Costa et al. (2000). The best fit physical pa- 
rameters of T b r = 2.5 x 10 4 K and n\V = 3.68 x 10 56 cur 3 
are obtained for the thermal bremsstrahlung emission. We 
note that according to the analysis of the IUE data to- 
wards pre-main-sequence stars by Johns-Krull, Valenti, & 
Linsky (2000), many classical T Tauri stars have color 
temperatures of ~ 10 4 K, derived from the mean con- 
tinuum flux at the wavelengths of 1958A and 1760A, in- 
dependent of the effective stellar temperatures. For the 
Ly a line emission, we adopt the central wavelength of 
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Fig. C.l. The radiation flux of the central T Tauri star 
(solid line) modeled by black body radiation (dashed line), 
thermal bremsstrahlung emission (dotted line), and Ly a 
line emission (dot-dashed line). The filled diamonds show 
the observational data towards TW Hya. 



Ao = 1215. 67A, the luminosity ratio of the line peak to the 
continuum radiation of LA ,Lya/(iA ,bb + ^A ,br) = 10 3 , 
and the line width of a = 2.0lA, following Herczeg et al. 
(2002; see also Ardila et al. 2002). 

The radiation flux from the central star is calculated 

as Fx, st ar = (£A,bb + £A,br + LA,Lya)/(47Ti? 2 ) and F A , st ar = 

iA,bb/(47ri? 2 ) for the models with and without UV excess 
radiation, respectively. We note that F^star in Sect. 2.2 is 
related to F\, s tar as AF^star = vFv, s t a T, where v is the fre- 
quency. In this paper we adopt the physical parameters of 
F*, T br , A , £A ,Lya/(iAo,bb + iA ,br), and a as mentioned 
above, and set the stellar radius to be 2i? Q for calculating 
the radiation from the central star, that is, we take the 
luminosity of the central star (L* = AiiR 2 F*) to be four 
times larger than that of TW Hya. 



Appendix D: Grain opacity 

In order to calculate the absorption (k„) and scattering 
(ov) coefficients of dust grains, we use the following dust 
model in this paper. First, as the dust components we 
adopt silicate and carboneous grains, both of which are 
considered as components of interstellar dust, and water 
ice, which is expected to be formed in the cold and dense 
region of protoplanetary disks. The optical properties of 
the carboneous grains are assumed to have a continuous 
distribution of graphite-like properties for larger sizes and 
PAH- like properties in the small size limit, following Li 
& Draine (2001). The mass fractional abundances of the 
above species are taken to be consistent with the solar 
elemental abundances: £ s u = 0.0043, Ccarbon = 0.0030, 
and (ice = 0.0094 (Anders & Grevesse 1989). Their bulk 
densities are set to be p S i\ — 3.5, /9 gr aphite = 2.24, and 
Pico = 0.92 g cm" 3 (see Li & Draine 2001 for the PAH's 
bulk density). Their sublimation temperatures are sim- 
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Fig. D.l. The monochromatic absorption coefficient 
(thick solid line) of dust grains consisting of carboneous 
(dashed line), silicate (thin solid line), and water ice (dot- 
dashed line). 

ply assumed to be T s n = 1500K, T ca rbon = 2300K, and 
T ico = 150K (e.g., Adams & Shu 1986). 

We assume that the silicate and carboneous grain par- 
ticles have the size distribution obtained by WDOla, which 
can reproduce the observational extinction curve of dense 
clouds with the ratio of visual extinction to reddening 
R v = A(V)/E(B - V) = 5.5 (see also Cardelli, Clayton 
& Mathis 1989). We use the model with the parameter of 
6c = 3.0 x 10~ 5 , which represents the total C abundance 
per H nucleus in the log-normal population, and the grain 
volumes are assumed to be the same as those of the dif- 
fuse interstellar medium (Case B of WDOla). The water 
ice is simply assumed to have the MRN size distribution of 
dn/da oc a -3 5 , where a is the radius of the dust particles 
(Mathis, Rumpl, & Nordsieck 1977). 

The total mass absorption coefficient k v and the total 
mass scattering coefficient u v are calculated as 

/ (dn/da) 8 a 3 {K„ t8 (a), a„ t8 (a)}da 

{ Ku ,a v } = Y, J - 7 . (D.l) 

s I (dn/da) s a 3 da 

where s represents the dust component - silicate, carbo- 
neous, or water ice. The mass absorption and scattering 
coefficients of dust particles with radius a is given by 

3 

{/«„ a (a), cr„ a (a)} = {Qabs(a, v), Q sca (a, v)}C, s (D.2) 

4ap s 

(e.g., Miyake & Nakagawa 1993). We adopt the absorp- 
tion factors Q a bs(a, v) described in Li & Draine (2001) for 
PAH molecules. The other absorption and scattering ef- 
ficiency factors, Qabs(a, v) and Q S ca{a,v), are computed 
by means of the Mie theory (Bohren & Huffman 1983), 
simply assuming spherical grains. In order to calculate 
the efficiency factors, we use the dielectric function of 
the 'smoothed UV astronomical silicate' and graphite by 
Draine & Lee (1984), Laor & Draine (1993), and WDOla 



(http://www.astro.prinston.edu/"" draine). Also we use 
the refractive indices by Miyake & Nakagawa (1993, see 
also references therein) for the water ice. 

The calculated monochromatic absorption coefficient 
is shown by the thick solid line in Fig. ID. II Each compo- 
nent of carboneous, silicate, and water ice is also plotted in 
the dashed, thin solid, and dot-dashed lines, respectively. 
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